from numpy.random import poisson
from vectors import *
import math
import random


def computeAttenuation(obs_Dict, total_dist):
	""" Computes the attentuation factor and the distance travelled through
	the obstacle medium. Returns the attenuation factor. Arguments are as 
	follows. 
	obstacles: Dictionary of Obstacles
	total_dist: Distance between Emitter and Detector"""
	if not (type(obs_Dict) is dict):
		raise TypeError("Expected a obstacle dictionary")
	
	attenuation = 0
	air_dist = float(total_dist)
	if not obs_Dict:
		for obs in obs_Dict:
			# Idx 1 of the obstacles contain distance of ray travelled
			air_dist -= obs_Dict[obs][1]
			# Changing cm^-1 to m^-1
			attenuation += Attenuations[obs.type][1] * \
							obs_Dict[obs][1] * 100
	
	# Add attenuation due to air
	attenuation += (Attenuations["Air"][1] * 100 * air_dist)

	#print "Attenuations: " + str(attenuation) + " " + str(air_dist)

	return attenuation

class PointEntity:
	"""The base class, from which GraphNode, Emitter, and Detector, and
	perhaps other point entities inherit from."""
	name = None
	def __init__(self, id, position, velocity=(0,0,0)):
		if not ((type(id) is int) or (type(id) is long)):
			raise TypeError("Expected integer or long for id, instead got " + \
			                str(type(id)))
		# it needs to be vec3-friendly or a pair of vec3-friendlies.
		if not (is_vec3_friendly(position) or \
		       (len(position) == 2 and \
		       is_vec3_friendly(position[0]) and \
		       is_vec3_friendly(position[1]))):
			raise ValueError("Expected something convertible to vec3 for " + \
			                 "position, instead got " + str(position))
		if not is_vec3_friendly(velocity):
			raise ValueError("Expected something convertible to vec3 for " + \
			                 "velocity, instead got " + str(velocity))
		self.id = id
		# if position is a box, indicating randomness,
		if len(position) == 2:
			corner1 = vec3(position[0])
			corner2 = vec3(position[1])
			self.p = corner1.min(corner2) + (abs(corner1 - corner2)).random()
		else:
			self.p = vec3(position)
		self.pos = self.p
		self.position = self.p
		self.v = vec3(velocity)
		self.vel = self.v
		self.velocity = self.v
		
	def __repr__(self):
		"""How the PointEntity sees itself."""
		# The idea here is to dig out the term PointEntity, or whatever other
		# inheriting class name it is, in order to only write __repr__ once.
		temp = str(self.__class__)
		namestr = ''
		if self.name:
			namestr = " (%s)" % self.name
		return ("%s %d" % (temp[temp.find('.')+1:], self.id)) + namestr
	
	def where(self, t=0.0):
		"""Returns the location of this entity, given a time relative to the
		present at which to compute position.  Assumes movement at the
		current velocity"""
		try:
			return self.p + self.v * float(t)
		except ValueError:
			raise ValueError(str(t) + " was not convertible to float.")

	def distanceTo(self, other):
		"""Returns the distance between this entity and another entity or
		vec3-like location"""
		if isinstance(other, PointEntity):
			return ((self.p - other.p) * (self.p - other.p))**.5
		if is_vec3_friendly(other):
			other_as_vec3 = vec3(other)
			return ((self.p - other_as_vec3) * (self.p - other_as_vec3))**.5
		raise ValueError("Expected a PointEntity or something " + \
		                 "convertible to vec3")


class Emitter(PointEntity):
	"""A point entity which emits non-negligible radiation of some kind."""
	### changed default to BozhilCs137 in the arguments, and intensity
	def __init__(self, id, position, velocity=(0,0,0), \
	             radType="BozhilCs137", intensity=1.0):
		try:
			self.intensity = float(intensity)
		except:
			raise TypeError("Expected float-friendly argument for intensity, " 
								+ "instead got a %s" % type(intensity))
		# Input Type Checking
		if not ((type(id) is int) or (type(id) is long)):
			raise TypeError("Expected integer or long for id, instead got " + \
			                str(type(id)))
		if not (is_vec3_friendly(position) or \
		       (len(position) == 2 and \
		       is_vec3_friendly(position[0]) and \
		       is_vec3_friendly(position[1]))):
			raise ValueError("Expected something convertible to vec3 for " + \
			                 "position, instead got " + str(position))
		if not is_vec3_friendly(velocity):
			raise ValueError("Expected something convertible to vec3 for " + \
			                 "velocity, instead got " + str(velocity))
		if not (type(radType) is str):
			raise TypeError("Expected a String as the Radiation type")

		self.id = id
		# if position is a box, indicating randomly variable location,
		if len(position) == 2:
			corner1 = vec3(position[0])
			corner2 = vec3(position[1])
			self.p = corner1.min(corner2) + (abs(corner1 - corner2)).random()
		else:
			self.p = vec3(position)
		self.pos = self.p
		self.position = self.p
		self.v = vec3(velocity)
		self.radiation = {}

		# References to the Simulation Object
		self.sim = None

		# Sets the type of the radiation when initializing
		self.setType(radType)
		### a quick access that might be useful long-term for other things
		self.quickCounts = sum(self.radiation["spectrum"]) * self.intensity

	def setIntensity(self, intensity):
		"""Sets the intensity to the new value, taking care to also set
		quickCounts to the appropriate value."""
		try:
			self.intensity = float(intensity)
		except:
			raise ValueError("Expected something float-friendly for an " + \
			                 "intensity value, instead got %s" % intensity)
		self.quickCounts = sum(self.radiation["spectrum"]) * self.intensity

	def setType(self, type):
		"""This function defines the emitter and its type. Must prove the type 
		as a string as the function input, and sets the radiation parameters
		as given by the type. Usually used in the __init__ function"""
		if type == "Cs137":
			self.radiation["type"] = "Cs137"
			self.radiation["num_energy"] = len(Cs137Energies)
			self.radiation["energy"] = Cs137Energies
			self.radiation["spectrum"] = Cs137Means[i]
		elif type == "BozhilCs137":
			self.radiation["type"] = "BozhilCs137"
			self.radiation["num_energy"] = len(Bozhil_Cs137_10keV)
			self.radiation["energy"] = Bozhil_Bins
			self.radiation["spectrum"] = Bozhil_Cs137_10keV
		else:
			raise ValueError("Requires the correct type of Radiation Source")


class Detector(PointEntity):
	"""A point entity which detects radiation.  Has an area and a dictionary of
	radiation bands it can detect, with detection intensities if necessary.
	That is, the entry "Cs":.75 in the dictionary would imply that Cesium-style
	radiation is detected at 75% of full intensity.
	
	NEW: The detector can be collimated so that it only receives radiation
	within a certain cone of the viewing direction.  The 'theta' attribute 
	refers to the angle between the direction the detector is facing and the
	surface of this cone.  (By default, theta is pi, i.e., the whole sphere)
	It is assumed to detect with perfect efficiency within the cone and to
	detect nothing outside the cone."""
	
	def __init__(self, id, position, velocity=(0,0,0), \
	             direction=(1,0,0)):
		if not ((type(id) is int) or (type(id) is long)):
			raise TypeError("Expected integer or long for id, instead got " + \
			                str(type(id)))
		if not is_vec3_friendly(position):
			raise ValueError("Expected something convertible to vec3 for " + \
			                 "position, instead got " + str(position))
		if not is_vec3_friendly(velocity):
			raise ValueError("Expected something convertible to vec3 for " + \
			                 "velocity, instead got " + str(velocity))
		if not is_vec3_friendly(direction):
			raise ValueError("Expected something convertible to vec3 for " + \
			                 "direction, instead got " + str(direction))
		self.theta = math.pi
		self.mode = "rad"
		self.id = id
		self.p = vec3(position)
		self.v = vec3(velocity)
		
		self.dir = vec3(direction).normalize()
		self.wedge = 0
		### Need to not hard-code this.
		self.num_wedges = 360
		
		wedge_arc = 2 * math.pi / self.num_wedges
		# north half
		if (self.dir.x >= 0 and self.dir.y >= 0) or \
		   (self.dir.x <  0 and self.dir.y > 0):
			angle = self.dir.angle(vec3(1,0,0))
			self.wedge = int(angle / wedge_arc)
		# south half
		else:
			angle = self.dir.angle(vec3(-1,0,0))
			self.wedge = int(angle / wedge_arc) + self.num_wedges/2
		
		self.radiation = {}
		self.hits = 0
		# the hits it received in the last timestep.
		self.recent_hits = 0

		# Detector Sensitivity, number of crystals of the detector
		self.size = 6
		
		# Simulation References
		self.sim = None

	def setTheta(self, theta=math.pi, mode="rad"):
		"""Sets the theta value of the detector, updating the arc fraction
		and window fraction as well."""
		try:
			self.theta = float(theta)
		except ValueError:
			raise ValueError("Expected something convertible to float for " + \
			                 "theta, instead got " + str(theta))
		if (mode != "rad") and (mode != "deg"):
			raise ValueError("Mode must be rad or deg (radians or degrees)")
		else:
			self.mode = mode

		if mode == "rad":
			self.arc_fraction = self.theta / math.pi
			self.window_fraction = (1 - math.cos(self.theta)) / 2
		else: # deg
			self.arc_fraction = self.theta / 180.0
			self.window_fraction = (1 - math.cos(self.theta * math.pi / 180)) / 2

#### Needs overhaul badly
	def setRotationalVelocity(self, value, mode="wedges"):
		"""Sets the rotational velocity of this detector.  Currently, this
		is only in wedges per timestep.  Which is terrible.
		TODO: Make this a nice, reasonable function."""
		if mode == "wedges":
			if type(value) is not int:
				print "Warning: casting value for setRotationalVelocity " + \
				      "as an integer, for wedges per timestep."
			self.wedges_per_timestep = int(value)
		else:
			print "that is not supported yet! defaulting to 1."
			self.wedges_per_timestep = 1

	def initializeWedgeDecays(self):
		"""Constructs a tuple of probability attenuations based on cosine
		of the angle from the center wedge.  That is to say, the first value
		is 1, the last is 0, and the ones in between decay as a
		quarter of the sinusoidal wave does.
		TODO: confirm method, write better docstring"""
		n = self.num_wedges / 4
		temp = [math.cos(i * math.pi / (2 * n)) for i in range(n)]
		self.wedge_decays = tuple(temp)

	def stepWedge(self):
		"""increments the wedge and direction."""
		self.wedge = (self.wedge + 1) % self.num_wedges
		self.dir = self.sim.heatmap.wedge2dir(self.wedge)

	def resetHits(self):
		"""Sets the hits variable, which is the count of photons that have
		struck this detector, back to zero."""
		self.hits = 0

	def getRadCounts(self, timestep, e, use_obs=True, use_cos=False, \
	                 bidirectional=False):
		"""Returns how many counts that it would get for the given time at the
		given input distance. Takes inputs as follows.
		time = time elapsed for calculating
		e = Emitter Object
		obstacles = dictionary of obstacles between obstacle and emitter."""
		if not (type(timestep) in (int, long, float)):
			raise TypeError("Requires a number as time argument")
		if not isinstance(e, Emitter):
			raise TypeError("Expected instance of Emitter for e, " + \
			                "instead got " + str(type(e)))
		if use_obs not in (False, True, 1, 0):
			raise TypeError("Requires a boolean as the optional argument")
		
		# Total Counts Made
		counts = 0
		dist = self.distanceTo(e)

		# Total attenuation factor
		# Calculating Obstacle Attenuation factors
		if use_obs:
			obstacles = self.sim.obstacleLocators[self]
			attenuation = computeAttenuation(obstacles, dist)
		else:
			attenuation = (Attenuations["Air"][1] * 100 * dist)
		
		if use_cos:
			incident = (e.p - self.p).normalize()
			# This is the value of the cosine of the angle to the normal or
			# antinormal, whichever is smaller.  They are effectively mirror
			# images (cosine is an even function)
			if bidirectional:
				reduction_due_to_angle = abs(incident * self.dir)
			else:
				reduction_due_to_angle = max(incident * self.dir, 0)
		else:
			reduction_due_to_angle = 1.0
		
		# Changed some calculation methods, all using data for 200keV

		mu = timestep * math.exp(-attenuation)/(dist ** 2)

		# For the given spectrum, return the total sum the poisson numbers
		# generated from the product of its element and mu	
		
		#for spec in e.radiation["spectrum"]:
		#	counts += poisson(spec * mu * self.size * e.intensity)
		
		# Initial Fix for Faster Simulation
		### revised by jeremy, note, quickCounts includes intensity
		counts = poisson(e.quickCounts * mu * self.size \
		                 * reduction_due_to_angle)

		self.recent_hits = counts
		return counts

	def getRadiation(self, time, e, obstacles=None):
		"""Returns the list of radiated particle energies for the given time
		and the given distance. TODO: Also need to adjust for obstacles!"""
		if not (type(time) in (int, long, float)):
			raise TypeError("Requires a number as time argument")
		if not isinstance(e, Emitter):
			raise TypeError("Expected instance of Emitter for e, " + \
			                "instead got " + str(type(e)))
		particles = []
		dist = self.distanceTo(e)

		if obstacles is not None:
			attenuation = computeAttenuation(obstacles, dist)
		else:
			attenuation = (Attenuations["Air"][1] * dist * 100)
		
		mu = time * math.exp(-dist * attenuation)/(dist ** 2)
		
		i = 0
		while i < e.radiation["num_energy"]:
			num = poisson(e.radiation["spectrum"][i] * mu * \
							 e.intensity * self.size)
			# If the Poisson number is <= 0, then it is not radiated particle
			if num <= 0:
				i += 1
				continue
			j = 0
			while j < num:
				particles.append(Bozhil_Bins[i] + random.random() * \
							(e.radiation["energy"][i] - Bozhil_Bins[i]))
				j += 1
			i += 1
			
		return particles

	def estimateIntensity(self, point):
		"""Given a location in space, the Detector estimates what intensity
		a source located there would have to be in order to produce the
		data it has received.
		mu = Lambda r^2 / (A * decay [* sweep fraction])
		Consider renaming, as there is an estimateIntensity in VoxelMap"""
		if is_vec3_friendly(point):
			p = vec3(point)
		elif isinstance(point, PointEntity):
			p = point.p
		else:
			raise ValueError("Expected a vec3-friendly or PointEntity as " + \
			                 "argument, instead got %s" % point)
		Lambda = float(self.hits) / (self.size * self.sim.time) \
		         - self.sim.assumed_background * self.window_fraction
		r =  self.distanceTo(p)
		# the A coefficient, assuming Bozhil data, is quickCounts / 1m air decay
		# i.e., .0159 == computeAttenuation({}, 1)
		A = 221.6394# * math.exp(-.0159)
		# this should be the product of exponential decays, i need a function
		# for this 
		obstacles = self.sim.obstaclemap.getObstaclesBetPts(self.p, p)
		attenuation = computeAttenuation(obstacles, r)
		decay = math.exp(-attenuation)
		# should arc_fraction come into play earlier?
		estimate = Lambda * r * r / (A * decay * self.arc_fraction)
		return estimate
		
		

class GraphNode(PointEntity):
	"""The foundation for the graph structure that will be used to represent
	the ability for entities to move throughout the system."""
	neighbors = 0

	def addNeighbor(self, graph, neighbor):
		"""Add another node to this node's set of neighbors. Must provide the
		graph dictionary as the function argument."""
		if type(graph) is not dict:
			raise TypeError("Require a graph dictionary as the input")
		if (neighbor > (len(graph) - 1)) or (neighbor < 0):
			raise ValueError("Cannot add a non-existing node")
		if neighbor in graph[self.id]:
			raise ValueError("Already a neighbor")
		else:
			# Check since no directed edges yet
			if self.id in graph[neighbor]:
				raise ValueError("Broken Graph!!! Debug")
			# Set neighbors on both sides
			graph[self.id] = graph[self.id] + (neighbor,)
			graph[neighbor] = graph[neighbor] + (self.id,)

	def removeNeighbor(self, graph, neighbor):
		"""Remove the selected node from this node's set of neighbors. 
		Must provide the graph dictionary as the function argument."""
		if type(graph) is not dict:
			raise TypeError("Require a graph dictionary as the input")
		if (neighbor > (len(graph) - 1)) or (neighbor < 0):
			raise ValueError("Cannot remove a non-existing node")
		if neighbor not in graph[self.id]:
			raise ValueError("Not a Neighbor of this node")
		else:
			# Check since no directed edges yet
			if self.id not in graph[neighbor]:
				raise ValueError("Broken Graph!!! Debug")
			# Remove neighbors from both sides
			temp = list(graph[self.id])
			temp.remove(neighbor)
			graph[self.id] = tuple(temp)
			temp = list(graph[neighbor])
			temp.remove(self.id)
			graph[neighbor] = tuple(temp)

	def isNeighbor(self, graph, neighbor):
		"""Checks to see if the input node is a neighbor of the current node.
		Must provide the graph dictionary as the function argument. """
		if type(graph) is not dict:
			raise TypeError("Require a graph dictionary as the input")
		if (neighbor > (len(graph) - 1)) or (neighbor < 0):
			raise ValueError("Cannot remove a non-existing node")
		if neighbor in graph[self.id]:
			# Check since no directed edges yet
			if self.id not in graph[neighbor]:
				raise ValueError("Broken Graph!!! Debug")
			else:
				return True
		
		return False


def isReachable(graph, id1, id2):
	'''Simple Depth-First search to see if there exists a path between 
	two nodes. Takes the adjacency matrix, and the ids of source and destination
	nodes.'''
	num_nodes = len(graph)
	if id1 < 0 or id2 < 0:
		raise ValueError("Input node id must be positive")
	if (id1 > (num_nodes - 1)) or (id2 > (num_nodes - 1)):
		raise ValueError("Node id exceeds the number of nodes")
	if id1 == id2:
		return True

	adj_Mat = adj_Mat.tolist()
	# Set of visited nodes
	visited = { id1: 1 }
	# Storing all neighbors
	neighbor_stack = []
	for i in range(num_nodes):
		if (adj_Mat[id1 - 1][i] == 1) and ((i + 1) != id1):
			neighbor_stack.append(i+1)

	while (len(neighbor_stack) > 0):
		neighbor = neighbor_stack.pop()
		# If the neighbor is our goal, we can reach him!
		if neighbor == id2:
			return True
		if neighbor not in visited:
			visited[neighbor] = 1
			for node in graph[neighbor]:
				neighbor_stack.append(node)
	return False

def shortestPath(graph, id1, id2):
	'''Uses Dijkstra's Shortest Path Algorithm to find the shortest path
	between two nodes. Takes the adjacency matrix, and the ids of source and
	destination nodes.'''
	num_nodes = len(graph)
	if id1 < 0 or id2 < 0:
		raise ValueError("Input node id must be positive")		
	if (id1 > (num_nodes - 1)) or (id2 > (num_nodes - 1)):		
		raise ValueError("Node id exceeds the number of nodes")		
		
	# Distances and Precedessors Dictionary		
	prev_dist = {}		
	final_dist = {}		
	prec = {}		
	# Sets all element to have no previous node and infinite distance		
	for node in graph:		
		prev_dist[node] = float('inf')		
		prec[node] = -1		
	# Source is 0 distance away		
	prev_dist[id1] = 0		
	final_dist[id1] = 0		
		
	while (len(prev_dist) > 0):		
		min_node = min(prev_dist, key=prev_dist.get)		
		if prev_dist[min_node] == float('inf'):		
			break		
		for neighbor in graph[min_node]:		
			temp_dist = prev_dist[min_node] + 1		
			if neighbor in final_dist:		
				if temp_dist < final_dist[neighbor]:		
					prev_dist[neighbor] = temp_dist		
					final_dist[neighbor] = temp_dist		
					prec[neighbor] = min_node		
			else:		
				prev_dist[neighbor] = temp_dist		
				final_dist[neighbor] = temp_dist		
				prec[neighbor] = min_node		
		del prev_dist[min_node]		

	#TODO: Return only the path when required	
	return (final_dist, prec)

# Linear Attenuation Coefficient(cm^(-1)) [100 keV, 200 keV, 500 keV]
# For Concrete, Density is assumed to be 2.3 g/cm^2 (usually in 2.2~2.4)
Attenuations = \
	{
		"Air": (0.000195, 0.000159, 0.000112),
		"Water": (0.167, 0.136, 0.097),
		"Carbon": (0.335, 0.274, 0.196),
		"Aluminum": (0.435, 0.324, 0.227),
		"Iron": (2.72, 1.09, 0.665),
		"Copper": (3.8, 1.309, 0.73),
		"Concrete": (0.39974, 0.29486, 0.205045)
	}


# keV energies (low edge bins) of the Cs137 spectrum
Cs137Energies = \
	(
		0.00,6.25,12.50,18.75,25.00,31.25,37.50,43.75,50.00,56.25,62.50,68.75,
		75.00,81.25,87.50,93.75,100.00,106.25,112.50,118.75,125.00,131.25,
		137.50,143.75,150.00,156.25,162.50,168.75,175.00,181.25,187.50,193.75,
		200.00,206.25,212.50,218.75,225.00,231.25,237.50,243.75,250.00,256.25,
		262.50,268.75,275.00,281.25,287.50,293.75,300.00,306.25,312.50,318.75,
		325.00,331.25,337.50,343.75,350.00,356.25,362.50,368.75,375.00,381.25,
		387.50,393.75,400.00,406.25,412.50,418.75,425.00,431.25,437.50,443.75,
		450.00,456.25,462.50,468.75,475.00,481.25,487.50,493.75,500.00,506.25,
		512.50,518.75,525.00,531.25,537.50,543.75,550.00,556.25,562.50,568.75,
		575.00,581.25,587.50,593.75,600.00,606.25,612.50,618.75,625.00,631.25,
		637.50,643.75,650.00,656.25,662.50,668.75,675.00,681.25,687.50,693.75,
		700.00,706.25,712.50,718.75,725.00,731.25,737.50,743.75,750.00,756.25,
		762.50,768.75,775.00,781.25,787.50,793.75,800.00 		
	)

	# Poisson mean number of photons counted per second for a detector at 1m 
	# distance from a 1 mCurie source
Cs137Means = \
	(
		1.8119,2.3141,8.2445,9.2376,13.8134,26.8366,9.1224,5.2549,5.0100,4.7809,
		4.707,4.6217,4.4559,4.3252,4.2409,4.2315,4.1257,4.2837,4.3881,4.3210,
		4.2748,4.1474,4.2189,4.2222,4.1937,4.1718,4.3327,4.2992,4.6091,5.3206,
		5.8525,4.8787,4.1808,3.7290,3.3584,3.2150,3.0018,2.8875,2.7341,2.6688,
		2.6358,2.4749,2.4042,2.3424,2.2581,2.2496,2.1426,2.0793,2.0236,2.0191,
		1.9244,1.8862,1.9142,1.7963,1.8150,1.7826,1.7773,1.7737,1.7309,1.7012,
		1.7327,1.7082,1.6804,1.6728,1.6804,1.6881,1.7192,1.7199,1.6693,1.6699,
		1.6395,1.6368,1.6155,1.6193,1.4460,1.3852,1.1026,0.7359,0.4895,0.4085,
		0.3765,0.3502,0.3201,0.3131,0.3195,0.2932,0.2827,0.3061,0.3050,0.2940,
		0.3061,0.2892,0.3253,0.3200,0.3360,0.3489,0.3844,0.3862,0.4440,0.5012,
		0.5487,0.6542,0.7793,0.9982,1.3826,2.5423,4.3900,1.8362,0.0315,0.0000,
		0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,
		0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000 		
	)
	
Bozhil_Cs137_10keV = (
		0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,
		0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,
		0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,
		0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0001,0.0000,
		0.0003,0.0003,0.0004,0.0007,0.0036,0.0050,0.0093,0.0161,0.0379,0.0707,0.1077,0.1791,0.2670,0.3271,0.3919,
		0.4335,0.4425,0.4637,0.4767,0.4838,0.4530,0.4043,0.3873,0.3733,0.3632,0.3828,0.3970,0.4034,0.4002,0.3930,
		0.4047,0.4227,0.4323,0.4356,0.4373,0.4583,0.4530,0.4764,0.4860,0.4692,0.4982,0.4894,0.5127,0.5062,0.5352,
		0.5392,0.5323,0.5267,0.5548,0.5630,0.5641,0.5609,0.5917,0.5884,0.5958,0.5893,0.5920,0.6260,0.6177,0.6286,
		0.6267,0.6427,0.6562,0.6491,0.6656,0.6853,0.6806,0.7098,0.7303,0.7205,0.7474,0.7307,0.7596,0.7901,0.7797,
		0.8025,0.8222,0.8276,0.8330,0.9018,0.8967,0.8989,0.9177,0.9498,0.9671,0.9961,1.0103,0.9978,1.0190,1.0144,
		1.0273,1.0460,1.0964,1.0898,1.0975,1.1090,1.1156,1.1082,1.1181,1.1047,1.1499,1.1546,1.1572,1.1594,1.1758,
		1.1766,1.2127,1.2268,1.2203,1.2265,1.2588,1.2830,1.2975,1.3493,1.3744,1.4380,1.5007,1.6212,1.7482,1.7883,
		1.8151,1.8744,1.8312,1.8094,1.7903,1.7955,1.7296,1.6756,1.6410,1.5953,1.5580,1.5029,1.4787,1.4170,1.4058,
		1.3949,1.3349,1.3012,1.2801,1.2544,1.2595,1.1858,1.1399,1.1280,1.1133,1.0839,1.0564,1.0360,1.0301,0.9820,
		0.9752,0.9432,0.9212,0.9254,0.9140,0.8923,0.8485,0.8453,0.8159,0.8012,0.7980,0.8011,0.7832,0.7583,0.7308,
		0.7258,0.7293,0.7139,0.7041,0.6815,0.6988,0.6666,0.6596,0.6728,0.6629,0.6562,0.6446,0.6576,0.6410,0.6470,
		0.6319,0.6278,0.6356,0.6296,0.6569,0.6094,0.6131,0.6247,0.6405,0.6161,0.6249,0.6406,0.6418,0.6231,0.6415,
		0.6206,0.6360,0.6403,0.6302,0.6445,0.6696,0.6603,0.6657,0.6841,0.6905,0.6854,0.7328,0.7177,0.7491,0.7459,
		0.7335,0.7876,0.7749,0.7920,0.8079,0.8087,0.8140,0.8214,0.8289,0.8394,0.8444,0.8113,0.7941,0.7687,0.6800,
		0.6391,0.5529,0.4739,0.4164,0.3638,0.3238,0.2835,0.2806,0.2462,0.2291,0.2128,0.2172,0.2107,0.1971,0.1824,
		0.1812,0.1728,0.1746,0.1630,0.1648,0.1614,0.1714,0.1637,0.1647,0.1659,0.1704,0.1746,0.1848,0.1824,0.1918,
		0.1945,0.1958,0.2142,0.2131,0.2326,0.2452,0.2668,0.2781,0.2939,0.3361,0.3805,0.4365,0.5125,0.6217,0.7477,
		0.9620,1.2304,1.6063,2.0473,2.6312,3.1686,3.5160,3.3308,2.7769,1.9974,1.3272,0.8297,0.5003,0.3107,0.1880,
		0.1239,0.0867,0.0648,0.0443,0.0356,0.0270,0.0194,0.0130,0.0089,0.0055,0.0041,0.0048,0.0027,0.0018,0.0011,
		0.0008,0.0002,0.0000,0.0000,0.0000,0.0007,0.0007,0.0000,0.0000,0.0000,0.0002,0.0000,0.0000,0.0000,0.0000,
		0.0000,0.0003,0.0000,0.0000,0.0000,0.0000,0.0000,0.0004,0.0002,0.0000,0.0000,0.0008,0.0000,0.0000,0.0000,
		0.0000,0.0000,0.0001,0.0000,0.0008,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,
		0.0000,0.0000,0.0011,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0002,
		0.0000,0.0000,0.0001,0.0000,0.0000,0.0000,0.0001,0.0000,0.0000,0.0000,0.0000,0.0007,0.0000,0.0000,0.0004,
		0.0000,0.0000,0.0004,0.0000,0.0000,0.0000,0.0008,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0003,0.0000,
		0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,
		0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0005,0.0000,0.0001,0.0000,0.0000,0.0000,0.0000,
		0.0001,0.0001,0.0000,0.0000,0.0000,0.0000,0.0000,0.0002,0.0000,0.0000,0.0000,0.0000,0.0002,0.0000,0.0000,
		0.0000,0.0000,0.0001,0.0001,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,
		0.0000,0.0000
	)
	
Bozhil_Background_10keV = (
		0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,
		0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,
		0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,
		0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,
		0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.0001,0.0002,0.0004,0.0007,0.0010,
		0.0019,0.0022,0.0031,0.0039,0.0046,0.0061,0.0060,0.0080,0.0075,0.0093,0.0100,0.0116,0.0115,0.0130,0.0144,
		0.0152,0.0161,0.0160,0.0181,0.0199,0.0213,0.0219,0.0240,0.0261,0.0276,0.0300,0.0317,0.0332,0.0338,0.0335,
		0.0354,0.0351,0.0361,0.0373,0.0386,0.0403,0.0416,0.0421,0.0436,0.0451,0.0471,0.0493,0.0483,0.0500,0.0497,
		0.0529,0.0530,0.0534,0.0567,0.0568,0.0569,0.0598,0.0597,0.0597,0.0606,0.0626,0.0633,0.0636,0.0648,0.0657,
		0.0670,0.0675,0.0668,0.0675,0.0696,0.0693,0.0676,0.0698,0.0712,0.0699,0.0707,0.0696,0.0735,0.0703,0.0708,
		0.0721,0.0700,0.0711,0.0705,0.0686,0.0688,0.0695,0.0680,0.0665,0.0671,0.0669,0.0672,0.0660,0.0645,0.0627,
		0.0613,0.0613,0.0609,0.0598,0.0597,0.0580,0.0571,0.0574,0.0557,0.0550,0.0538,0.0536,0.0533,0.0517,0.0509,
		0.0503,0.0479,0.0484,0.0487,0.0468,0.0463,0.0456,0.0440,0.0436,0.0430,0.0433,0.0408,0.0409,0.0401,0.0384,
		0.0384,0.0377,0.0366,0.0379,0.0361,0.0361,0.0359,0.0363,0.0359,0.0342,0.0337,0.0333,0.0311,0.0288,0.0286,
		0.0281,0.0269,0.0263,0.0255,0.0249,0.0243,0.0233,0.0235,0.0226,0.0218,0.0223,0.0219,0.0208,0.0198,0.0205,
		0.0207,0.0198,0.0197,0.0200,0.0189,0.0184,0.0183,0.0179,0.0164,0.0168,0.0157,0.0163,0.0160,0.0155,0.0151,
		0.0148,0.0147,0.0148,0.0148,0.0143,0.0145,0.0149,0.0146,0.0146,0.0143,0.0144,0.0130,0.0119,0.0119,0.0120,
		0.0109,0.0111,0.0111,0.0108,0.0099,0.0102,0.0099,0.0100,0.0093,0.0098,0.0095,0.0087,0.0094,0.0086,0.0085,
		0.0093,0.0083,0.0083,0.0087,0.0075,0.0083,0.0087,0.0079,0.0081,0.0078,0.0081,0.0076,0.0077,0.0075,0.0074,
		0.0070,0.0069,0.0066,0.0068,0.0066,0.0065,0.0068,0.0068,0.0062,0.0065,0.0063,0.0065,0.0070,0.0068,0.0066,
		0.0067,0.0061,0.0061,0.0059,0.0055,0.0055,0.0061,0.0057,0.0052,0.0058,0.0051,0.0051,0.0055,0.0054,0.0055,
		0.0059,0.0055,0.0058,0.0054,0.0054,0.0053,0.0058,0.0055,0.0061,0.0057,0.0061,0.0054,0.0056,0.0048,0.0048,
		0.0047,0.0047,0.0044,0.0042,0.0045,0.0048,0.0043,0.0048,0.0050,0.0048,0.0048,0.0041,0.0044,0.0040,0.0041,
		0.0039,0.0036,0.0038,0.0036,0.0037,0.0037,0.0034,0.0035,0.0036,0.0038,0.0035,0.0038,0.0036,0.0035,0.0036,
		0.0034,0.0034,0.0029,0.0033,0.0034,0.0029,0.0029,0.0035,0.0033,0.0031,0.0028,0.0034,0.0032,0.0029,0.0031,
		0.0028,0.0028,0.0031,0.0031,0.0028,0.0030,0.0027,0.0028,0.0027,0.0029,0.0029,0.0026,0.0027,0.0028,0.0027,
		0.0027,0.0027,0.0028,0.0029,0.0024,0.0028,0.0031,0.0030,0.0028,0.0026,0.0026,0.0028,0.0027,0.0027,0.0029,
		0.0026,0.0023,0.0025,0.0030,0.0025,0.0028,0.0025,0.0024,0.0024,0.0025,0.0026,0.0026,0.0022,0.0022,0.0026,
		0.0023,0.0023,0.0023,0.0023,0.0023,0.0020,0.0023,0.0023,0.0020,0.0021,0.0024,0.0019,0.0021,0.0019,0.0021,
		0.0023,0.0023,0.0021,0.0023,0.0024,0.0024,0.0022,0.0021,0.0024,0.0023,0.0023,0.0019,0.0021,0.0022,0.0025,
		0.0025,0.0021,0.0022,0.0019,0.0020,0.0019,0.0021,0.0019,0.0019,0.0023,0.0021,0.0018,0.0016,0.0017,0.0017,
		0.0014,0.0014,0.0012,0.0009,0.0011,0.0011,0.0010,0.0011,0.0009,0.0011,0.0010,0.0010,0.0010,0.0009,0.0010,
		0.0010,0.0007,0.0009,0.0011,0.0008,0.0008,0.0009,0.0009,0.0008,0.0008,0.0008,0.0014,0.0010,0.0010,0.0010,
		0.0012,0.0013,0.0013,0.0016,0.0017,0.0021,0.0022,0.0022,0.0022,0.0022,0.0017,0.0014,0.0011,0.0008,0.0008,
		0.0009,0.0007
	)
	
Bozhil_Bins = (
		0.0057,0.0229,0.0515,0.0916,0.1431,0.2060,0.2804,0.3662,0.4635,0.5722,0.6924,0.8240,0.9670,1.1215,1.2875,
		1.4648,1.6537,1.8539,2.0657,2.2888,2.5234,2.7695,3.0270,3.2959,3.5763,3.8681,4.1714,4.4861,4.8122,5.1498,
		5.4989,5.8594,6.2313,6.6147,7.0095,7.4158,7.8335,8.2626,8.7032,9.1553,9.6188,10.0937,10.5801,11.0779,
		11.5871,12.1078,12.6400,13.1836,13.7386,14.3051,14.8830,15.4724,16.0732,16.6855,17.3092,17.9443,18.5909,
		19.2490,19.9184,20.5994,21.2917,21.9955,22.7108,23.4375,24.1756,24.9252,25.6863,26.4587,27.2427,28.0380,
		28.8448,29.6631,30.4928,31.3339,32.1865,33.0505,33.9260,34.8129,35.7113,36.6211,37.5423,38.4750,39.4192,
		40.3748,41.3418,42.3203,43.3102,44.3115,45.3243,46.3486,47.3843,48.4314,49.4900,50.5600,51.6415,52.7344,
		53.8387,54.9545,56.0818,57.2205,58.3706,59.5322,60.7052,61.8896,63.0856,64.2929,65.5117,66.7419,67.9836,
		69.2368,70.5013,71.7773,73.0648,74.3637,75.6741,76.9958,78.3291,79.6738,81.0299,82.3975,83.7765,85.1669,
		86.5688,87.9822,89.4070,90.8432,92.2909,93.7500,95.2206,96.7026,98.1960,99.7009,101.2173,102.7451,104.2843,
		105.8350,107.3971,108.9706,110.5556,112.1521,113.7600,115.3793,117.0101,118.6523,120.3060,121.9711,123.6477,
		125.3357,127.0351,128.7460,130.4684,132.2021,133.9474,135.7040,137.4722,139.2517,141.0427,142.8452,144.6590,
		146.4844,148.3212,150.1694,152.0290,153.9001,155.7827,157.6767,159.5821,161.4990,163.4274,165.3671,167.3183,
		169.2810,171.2551,173.2407,175.2377,177.2461,179.2660,181.2973,183.3401,185.3943,187.4599,189.5370,191.6256,
		193.7256,195.8370,197.9599,200.0942,202.2400,204.3972,206.5659,208.7460,210.9375,213.1405,215.3549,217.5808,
		219.8181,222.0669,224.3271,226.5987,228.8818,231.1764,233.4824,235.7998,238.1287,240.4690,242.8207,245.1839,
		247.5586,249.9447,252.3422,254.7512,257.1716,259.6035,262.0468,264.5016,266.9678,269.4454,271.9345,274.4350,
		276.9470,279.4704,282.0053,284.5516,287.1094,289.6786,292.2592,294.8513,297.4548,300.0698,302.6962,305.3341,
		307.9834,310.6441,313.3163,316.0000,318.6951,321.4016,324.1196,326.8490,329.5898,332.3421,335.1059,337.8811,
		340.6677,343.4658,346.2753,349.0963,351.9287,354.7726,357.6279,360.4946,363.3728,366.2624,369.1635,372.0760,
		375.0000,377.9354,380.8823,383.8406,386.8103,389.7915,392.7841,395.7882,398.8037,401.8307,404.8691,407.9189,
		410.9802,414.0530,417.1371,420.2328,423.3398,426.4584,429.5883,432.7297,435.8826,439.0469,442.2226,445.4098,
		448.6084,451.8185,455.0400,458.2729,461.5173,464.7732,468.0405,471.3192,474.6094,477.9110,481.2241,484.5486,
		487.8845,491.2319,494.5908,497.9610,501.3428,504.7359,508.1406,511.5566,514.9841,518.4231,521.8735,525.3353,
		528.8086,532.2933,535.7895,539.2971,542.8162,546.3467,549.8886,553.4420,557.0068,560.5831,564.1708,567.7700,
		571.3806,575.0027,578.6362,582.2811,585.9375,589.6053,593.2846,596.9753,600.6775,604.3911,608.1161,611.8526,
		615.6006,619.3600,623.1308,626.9131,630.7068,634.5119,638.3286,642.1566,645.9961,649.8470,653.7094,657.5832,
		661.4685,665.3652,669.2734,673.1930,677.1240,681.0665,685.0204,688.9858,692.9626,696.9509,700.9506,704.9618,
		708.9844,713.0184,717.0639,721.1208,725.1892,729.2690,733.3603,737.4630,741.5771,745.7027,749.8398,753.9883,
		758.1482,762.3196,766.5024,770.6966,774.9023,779.1195,783.3481,787.5881,791.8396,796.1025,800.3769,804.6627,
		808.9600,813.2687,817.5888,821.9204,826.2634,830.6179,834.9838,839.3612,843.7500,848.1503,852.5620,856.9851,
		861.4197,865.8657,870.3232,874.7921,879.2725,883.7643,888.2675,892.7822,897.3083,901.8459,906.3950,910.9554,
		915.5273,920.1107,924.7055,929.3118,933.9294,938.5586,943.1992,947.8512,952.5146,957.1896,961.8759,966.5737,
		971.2830,976.0036,980.7358,985.4794,990.2344,995.0008,999.7787,1004.5681,1009.3689,1014.1811,1019.0048,
		1023.8400,1028.6865,1033.5445,1038.4140,1043.2949,1048.1873,1053.0910,1058.0063,1062.9330,1067.8711,1072.8207,
		1077.7817,1082.7541,1087.7380,1092.7334,1097.7402,1102.7584,1107.7881,1112.8292,1117.8818,1122.9458,1128.0212,
		1133.1081,1138.2065,1143.3163,1148.4375,1153.5702,1158.7143,1163.8699,1169.0369,1174.2153,1179.4052,1184.6066,
		1189.8193,1195.0436,1200.2792,1205.5264,1210.7849,1216.0549,1221.3364,1226.6293,1231.9336,1237.2494,1242.5766,
		1247.9153,1253.2654,1258.6269,1263.9999,1269.3844,1274.7803,1280.1876,1285.6064,1291.0366,1296.4783,1301.9314,
		1307.3959,1312.8719,1318.3594,1323.8583,1329.3686,1334.8904,1340.4236,1345.9682,1351.5244,1357.0919,1362.6709,
		1368.2613,1373.8632,1379.4765,1385.1013,1390.7375,1396.3852,1402.0443,1407.7148,1413.3968,1419.0903,1424.7952,
		1430.5115,1436.2392,1441.9785,1447.7291,1453.4912,1459.2648,1465.0497,1470.8462,1476.6541,1482.4734,1488.3041,
		1494.1463,1500.0000
	)
	
